ohibc logo
OHI British Columbia | OHI Science | Citation policy

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
library(data.table)
library(sf)


source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')  ###
  ### includes library(tidyverse); library(stringr); dir_M points to ohi directory

dir_git <- '~/github/spp_risk_dists'

### goal specific folders and info
dir_setup <- file.path(dir_git, 'data_setup')
dir_data  <- file.path(dir_git, 'data')
dir_o_anx <- file.path(dir_O, 'git-annex/spp_risk_dists')

source(file.path(dir_git, 'data_setup/api_fxns.R'))

### Gall-Peters doesn't have an EPSG?
gp_proj4 <- '+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'

1 Summary

Create a set of maps of the distribution of biodiversity intactness - all species assessed and mapped by IUCN. These maps are generated at 10 km2 resolution in a Gall-Peters projection. These maps will be generated using all available species:

  • Mean risk
  • Variance of risk
  • Number of species for mean/var calculations
  • Number of species categorized as “threatened” (i.e. VU, EN, CR)
  • Mean trend
  • Number of species for trend calculations

A selection of these maps will be generated for taxonomic groups and range sizes in a separate Rmd.

Future iterations may include:

  • Range-rarity-weighted mean and variance of risk
  • Range rarity-weighted species richness

2 Data Sources

IUCN Red List spatial data download IUCN Red List API Gina Ralph (IUCN)

3 Methods

3.1 Spatial distribution of current extinction risk

From the 1a and 1b biodiversity risk map scripts, gather the rasters for the various maps:

  • Mean risk (unweighted and range-rarity-weighted)
  • Pct threatened (unweighted and range-rarity-weighted)
agg_fact <- 1 ### for map creation, aggregate for speed and ease; but
              ### not for density plots!

risk_un_rast    <- raster(file.path(dir_git, 'output', 'mean_risk_raster_gp.tif'))
risk_rr_rast    <- raster(file.path(dir_git, 'output', 'mean_rr_risk_raster_gp.tif'))
threat_un_rast  <- raster(file.path(dir_git, 'output', 'pct_threat_raster_gp.tif'))
threat_rr_rast  <- raster(file.path(dir_git, 'output', 'sr_rr_pct_threat_raster_gp.tif'))
n_spp_rast <- raster(file.path(dir_git, 'output', 'n_spp_risk_raster_gp.tif'))
n_rr_rast  <- raster(file.path(dir_git, 'output', 'sr_rr_risk_raster_gp.tif'))

3.1.1 And now, the maps

First do a facet grid of mean risk (unweighted and rr-weighted); then do a facet grid of density; then combine with gridExtra or cowplot

four_panel <- function(map1_rast, map2_rast,
                       colors = c('green4', 'lightyellow', 'orange', 'red2', 'red4', 'purple4'),
                       values = c(0, .2, .4, .6, .8, 1.0),
                       limits = c(0, 1),
                       labels = c('LC', 'NT', 'VU', 'EN', 'CR', 'EX'),
                       breaks = c( 0.0,  0.2,  0.4,  0.6,  0.8,  1.0)) {

  land_poly <- sf::read_sf(file.path(dir_git, 'spatial/ne_10m_land/ne_10m_land.shp')) %>%
    st_transform(gp_proj4)
  
  map_df <- data.frame(val_1 = values(map1_rast),
                       val_2 = values(map2_rast)) %>%
    cbind(coordinates(map1_rast)) %>% 
    filter(!is.na(val_1))
  
  marg_val <- .25 ### plot.margin values in cm, with 0 between the map and density
  
  map1 <- ggplot(map_df) +
    geom_raster(aes(x, y, fill = val_2), show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map() +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = colors, values = values, limits = limits,
                         labels = labels, breaks = breaks,
                         guide  = guide_colourbar(label.position = 'left',
                                                  label.hjust = 1)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0))
  
  # ggsave('map1_test.png')
  
  map2 <- ggplot(map_df) +
    geom_raster(aes(x, y, fill = val_2), show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map() +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = colors, values = values, limits = limits,
                         labels = labels, breaks = breaks,
                         guide  = guide_colourbar(label.position = 'left',
                                                  label.hjust = 1)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0))

  ###################
  
  ### set up a dataframe of values to craft a color bar using geom_segment
  colorbar_df <- data.frame(x = seq(0, 1, .001), y = -1)

  dens1 <- ggplot(map_df) +
    ggtheme_plot() +
    geom_vline(xintercept = mean(map_df$val_1, na.rm = TRUE)) +
    geom_segment(data = colorbar_df, 
                 aes(x = x, xend = x, color = x), 
                 y = 0, yend = -10, size = 1,
                 show.legend = FALSE) +
    geom_density(aes(x = val_1, ..scaled..), alpha = .3, size = .25, fill = 'grey30') +
    scale_color_gradientn(colors = colors, values = values, limits = limits,
                          labels = labels, breaks = breaks) +
    theme(axis.text.x  = element_blank(), 
          axis.title = element_blank(),
          panel.grid.major.x = element_blank()) +
    scale_x_continuous(labels = labels, breaks = breaks, limits = limits,
                       expand = c(0, 0)) +
    scale_y_continuous(limits = c(-.15, 1)) +
    coord_flip()
    

  dens2 <- ggplot(map_df, aes(x = val_2)) +
    ggtheme_plot() +
    geom_vline(xintercept = mean(map_df$val_2, na.rm = TRUE)) +
    geom_segment(data = colorbar_df, 
                 aes(x = x, xend = x, color = x), 
                 y = 0, yend = -10, size = 1,
                 show.legend = FALSE) +
    geom_density(aes(x = val_2, ..scaled..), alpha = .3, size = .25, fill = 'grey30') +
    scale_color_gradientn(colors = colors, values = values, limits = limits,
                          labels = labels, breaks = breaks) +
    theme(axis.text.x  = element_blank(), 
          axis.title = element_blank(),
          panel.grid.major.x = element_blank()) +
    scale_x_continuous(labels = labels, breaks = breaks, limits = limits,
                       expand = c(0, 0)) +
    scale_y_continuous(limits = c(-.15, 1)) +
    coord_flip()
 
  panel_top <- cowplot::plot_grid(map1, dens1, 
                                  axis = 'b',
                                  rel_widths = c(4, 1))
  panel_btm <- cowplot::plot_grid(map2, dens2, 
                                  axis = 'b',
                                  rel_widths = c(4, 1))
  
  four_panel <- cowplot::plot_grid(panel_top, panel_btm, 
                                   labels = c('A', 'B'),
                                   nrow = 2, ncol = 1, align = 'v')

}
mean_four_panel <- four_panel(risk_un_rast, risk_rr_rast)
  ### function defaults set up for this plot
fig1_path <- file.path(dir_git, 'ms_figures/fig1_mean_maps_with_dens.png')
ggsave(plot = mean_four_panel, 
       filename = fig1_path, 
       width = 6, height = 6, units = 'in', dpi = 300)


threat_four_panel <- four_panel(threat_un_rast, threat_rr_rast,
                                colors = c('grey90', 'coral', 'red2', 'red4', 'purple4'),
                                values = c(0, .3, .5, .7, 1),
                                breaks = c(0, .25, .5, .75, 1),
                                labels = c('0%', '25%', '50%', '75%', '100%'))
fig2_path <- file.path(dir_git, 'ms_figures/fig2_pct_thr_maps_with_dens.png')
ggsave(plot = threat_four_panel, 
       filename = fig2_path,
       width = 6, height = 6, dpi = 300)

3.1.2 Fig 1: Biodiversity risk

3.1.3 Fig 2: Percent of species threatened

```